Qualitative and quantitative analysis of stability and instability dynamics of positive lattice solitons 
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We present a unified approach for qualitative and quantitative analysis of stability and instability dynamics 
of positive bright solitons in multi-dimensional focusing nonlinear media with a potential (lattice), which can 
be periodic, periodic with defects, quasiperiodic, single waveguide, etc. We show that when the soliton is 
unstable, the type of instability dynamic that develops depends on which of two stability conditions is violated. 
Specifically, violation of the slope condition leads to a focusing instability, whereas violation of the spectral 
condition leads to a drift instability. We also present a quantitative approach that allows to predict the stability 
and instability strength. 

PACS numbers: 42.65 Jx, 42.65 Tg, 03.75 Lm 



I. INTRODUCTION 

Solitons, or solitary waves, are localized nonlinear waves 
that maintain their shape during propagation. They are preva- 
lent in many branches of physics, and their properties have 
provided deep insight into complex nonlinear systems. The 
stability properties of solitons are of fundamental importance. 
Stable solitons are both natural carriers of energy in naturally 
occurring systems and often the preferred carriers of energy 
in engineered systems. Their stability also makes them most 
accessible to experimental observation. 

The first studies considered stability of solitons in homo- 
geneous media. In recent years there has been a considerable 
interest in the study of solitons in lattice-type systems. Such 
solitons have been observed in optics using waveguide arrays, 
photo-refractive materials, photonic crystal fibers, etc., in both 
one-dimensional and multidimensional lattices, mostly peri- 
odic sinusoidal square lattices C D, B [E B IH IE S Htl or sin- 
gle waveguide potentials |9l UOl [Till , but also in discontinu- 
ous lattices (surface solitons) 1121 . radially-symmetric Bessel 
lattices IU3I1 . lattices with triangular or hexagonal symme- 
try IBS lattices with defectsjQ^ynflMljMEEIll, 
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with quasicrystal structures lfl6ll23 , 24 , 25, 
random potentials ll29[|30ll . Solitons have also been observed 
in the context of Bose-Einstein Condensates (BEC) j3ll[32Tl . 
where lattices have been induced using a variety of tech- 
niques. 

Stability of lattice solitons has been studied in hundreds of 
papers. The majority of these papers focused on one specific 
physical configuration, i.e., a specific dimension (mostly in 
ID), nonlinearity and lattice type. In addition, in several stud- 
ies, general conditions for stability and instability were de- 
rived (see Section|III|. In all of these studies, the key question 
was whether the soliton is stable (yes) or unstable (no). 

Fibich, Sivan and Weinstein went beyond this binary view 
by developing a qualitative and quantitative approach to 
stability of positive lattice solitons. This was first carried 
out for spatially non-homogeneous nonlinear potentials in 
11331 13411 . These ideas were then developed by Sivan, Fibich 
and coworkers in the context of linear non-homogeneous po- 



tentials in 1 35, 361I37D- These studies showed that the quali- 
tative nature of the instability dynamics is determined by the 
particular violated stability condition. In addition, they pre- 
sented a quantitative approach for prediction of the stability 
or instability strength. Specifically, these papers considered 
the cases of a one-dimensional nonlinear lattice lO . a two- 
dimensional nonlinear lattice B34I1 . a one-dimensional linear 
delta-function potential ll36ll and narrow solitons in a linear 
lattice JH]. 

In the present article, the results of JH [M HI IH [H are 
combined into a unified theory for stability and instability of 
lattice solitons that can be summarized in a few rules (Sec- 
tion PVU l. We illustrate how these rules can be applied in a 
variety of examples that may be useful to experimental stud- 
ies. 



II. MODEL, NOTATION AND DEFINITIONS 

We study the stability and instability dynamics of lattice 
solitons of the nonlinear Schrodinger (NLS) equation with an 
external potential, which in dimensionless form is given by 

iA z (x, z) + AA + (1 - V nl (x)) F (\A\ 2 ) A - V t (x)A = 0. 

(1) 

Equation (fl} is also referred to as the Gross-Pitaevskii equa- 
tion (GP). NLS/GP underlies many models of nonlinear wave 
propagation in nonlinear optics and macroscopic quantum 
systems (BEC). For example, in the context of laser beam 
propagation, A(x, z) corresponds to the electric field ampli- 
tude, z > is the distance along the direction of propaga- 
tion, x = (xi, . . . , Xd) is the transverse d-dimensional space 
[e.g., the (x, y) plane for propagating in bulk medium] and 



A 



-az is the d-dimensional diffraction term. The 



nonlinear term models the intensity-dependence of the refrac- 
tive index. For example, = \A\ 2 corresponds to the 
optical Kerr effect and F {\A\ 2 ) = 1/ (l + |.A| 2 ) corresponds 
to photorefractive materials, see e.g., 13811 . The potentials V; 
and V n i correspond to a modulation of the linear and nonlinear 
refractive indices, respectively. In BEC, z = t is time, A(x, t) 
represents the wave function of the mean-field atomic conden- 
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sate, F (|A| 2 ) = \A\ 2 represents contact (cubic) interaction, 
and the potentials Vi(x) and V n i(x) are induced by externally 
applied electro-magnetic fields IB9I1 . 

We define a soliton to be any solution of Eq. (Q} of the form 
A(x, z) — u(x)e~' lfJ - z , where fi is the propagation constant 
and u(x), the soliton profile, is a real-valued function that de- 
cays to zero at infinity and satisfies 

Au+ (1- V n i(x)) F{u 2 )u + fiu - V t u = 0. (2) 

Solitons can exist only for /i in the gaps in the spectrum of 
the linear problem 

Au + (m — Viu = 0, (3) 

i.e., for values of fi such that the linear problem (0 does not 
have any non-trivial solution, see e.g., 114011 . 

Solitons in a lattice potential, or more general non- 
homogeneous potential, may be understood as bounds states 
of an effective (self-consistent) potential, V e ff = Vi(x) + 
(— 1 + V n i(x)) F (u 2 (x)). They arise (i) via bifurcation from 
the zero-amplitude state with energy at an end point of a 
continuous spectral band (finite or semi-infinite) of extended 
states of the linear operator of — A + Vi [80] or (ii) if Vi is a 
potential with a defect, via bifurcation from discrete eigenval- 
ues (localized linear modes) within the spectral gaps (semi- 
infinite or finite), which in addition to the bare nonlinearity, 
can serve to nucleate a localized nonlinear bound state 14711 . 

In this paper, we only consider positive solitons (u > 0) 
of both type (i) and (ii). This is always the case for the least 
energy state within the semi-infinite gap, i.e., when — oo < 
fi < /-i^2, where /i^2 is the lowest point in the spectrum of 
Eq. (0, at which the first band begins. Solitons whose fre- 
quencies lie in finite spectral gaps are usually referred to as 
gap solitons. However, gap solitons typically oscillate and 
change sign and are therefore not covered by the theory pre- 
sented in this paper [ 8811 . 

We study the dynamics of NLS/GP and its solitons in the 
space H 1 , with norm \\f\\ 2 m := J(|/| 2 + | V/| 2 |)df. The 
natural notion of stability is orbital stability, defined as fol- 
lows: 

Definition II.l Let u(x) be a solution of Eq. (0 with prop- 
agation constant fi. Then, the soliton solution u(x)e~ lfJ,z 
of NLS/GP eqn. (0 is orbitally stable if for all £ > 0, 
there exists 5(e) > such that for any initial condition 
Ao with inf 7G R \\Aq — iie* 7 ||^fi < 5, then for all 
z > the corresponding solution A(x, z) of Eq. (0 satisfies 
inf 7GR \\A(-, z) - ue^Wui < e . 

In discussing the stability theory for NLS it is useful to refer 
to its Hamiltonian structure: idtA = 5H/SA*, where 

H[A,A*} = J (\VA\ 2 + V t \A\ 2 -{l-V nl )G{\A\ 2 ))dx, 

and G'{s) = F(s), G(0) = 0. The Hamiltonian, H, and the 
optical power (particle number): 



are conserved integrals for NLS. Eq. (O for u(x; /i), the soli- 
ton profile, can be written equivalently as the energy station- 
arity condition 5£ = 0, where £ = H — [iP. Soliton stability 
requires a study of S 2 £, the second variational derivative of £ 
about u. For NLS in general, stable solitons need to be local 
energy minimizers; see e.g. M~k and also ll89ll . 



III. SOLITON STABILITY - OVERVIEW 

The first analytic result on soliton stability was obtained by 
Vakhitov and Kolokolov JU]. They proved, via a study of the 
linearized perturbation equation, that a necessary condition 
for stability of the soliton u(x; fi) is 

dP(a) 

—T 1 < 0, (4) 
a/i 

i.e. the soliton is stable only if its power decreases with 
increasing propagation constant /i. This condition will be 
henceforth called the slope condition. 

Subsequent studies of nonlinear stability analysis of soli- 
tons revealed the central role played by the number of negative 
and zero eigenvalues of the operator 

L+ = -A + Vi-{l- Vm) (F{u 2 ) - 2u 2 F'{u)) - (x, (5) 

which is the real part of S 2 £ J90ll . Weinstein [43] showed that 
for a homogeneous (translation-invariant) medium, (Vi = 0, 
V n i = 0), for u(x; n) > 0, if the slope condition © is satis- 
fied and L + has only one negative eigenvalue, then the soli- 
tons are nonlinearly stable. Later, in 14711 (Theorem 3.1; see 
also Theorem 6 of B48I1 ) it was shown that in the presence 
of a linear potential which is bounded below and decaying 
at infinity, solitons are stable if in addition, L + also has no 
zero eigenvalue(s). A related treatment was given to the nar- 
row soliton (semi-classical limit) subcritical nonlinearity case 
in i35l l44l |49[] and for solitons in spatially varying nonlinear 
potentials in 11331 13411 . 

General sufficient conditions for instability were given by 
Grillakis 15011 and Jones [51]. These results imply that if ei- 
ther the slope is positive or if L + has more than one negative 
eigenvalue, then the soliton is unstable. 

A direct consequence of the arguments of Section 3 and 
Theorem 3.1 in H, B, and iM El, is a stability theo- 
rem (used in this paper), which applies to positive solitons of 
NLS ([1]), whose frequencies lie in the semi-infinite spectral 
gap of —A + V; see also l53ll . 

Theorem III.l Let u(x) be a positive solution ofEq. (O with 
propagation constant fi within the semi-infinite gap, i.e. fi < 
I^be- Then, A = u(x)e~ tliz is an orbitally-stable solution of 
the NLS (Q if both of the following conditions hold: 

1. The sfope (Vakhitov -Kolokolov) condition: 



P = / \A\ 2 dx 



(6) 
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2. The spectral condition: L + has no zero eigenvalues 
and 

n_(L+) = l. (7) 

If either j£ > or ri- > 2, f/ze soliton is unstable. 
We note that Theorem lHI. ll does not cover two cases: 

1 . ^ = 0: For homogeneous media, V/ = V n i = 0, soli- 
tons are unstable; see B55l 16111 for power nonlinearities 
and lf52ll for general nonlinearities. There are no ana- 
lytic results for inhomogeneous media. 

2. n_(£+) = 1 and zero is an eigenvalue of multiplic- 
ity one or higher: This case will be discussed in Sec- 
tions|mA]and|IirB] 

We also note that Grillakis, Shatah and Strauss (GSS) H 
l45ll gave an alternative abstract formulation of a stability the- 
ory for positive solitons Hamiltonian systems, including NLS 
with a general class of linear and nonlinear spatially depen- 
dent potentials. In this formulation, the spectral condition on 
n_ (£+) and the slope condition are coupled, see detailed dis- 
cussion in [35]. The formulation of Theorem lHI. 11 is a more 
refined and stronger statement. Specifically, it decouples the 
slope condition and the spectral condition on n_ (£+) as two 
independent necessary conditions for stability and shows that 
a violation of either of them would lead to instability. This 
decoupling is at the heart of our qualitative approach since vi- 
olation of each condition leads to a different type of instability. 
Stability of solitons in homogeneous media has also been in- 
vestigated using the Hamiltonian-Power curves, see e.g., M46I1 . 

A. Review of stability conditions in homogeneous media 

Stability and instability of solitons in homogeneous media 
(i.e., V = 0) have been extensively investigated llsill . In this 
case, t^BE 

: ' = 0, i.e., the semi-infinite gap associated with 
Eq. (01 is (—oo,0). For every /j < and Xq S R d , there 
exists a soliton centered at Xq which is radially-symmetric in 
r = \x — xq\, positive, and monotonically decaying in r. 

In the case of a power-law nonlinearity F(|u|) = \u\ 2cr , the 
slope condition (O depends on the dimension d and nonlin- 
earity exponent a as follows 043115511 : 

1. In the subcritical case d < 2/cr, ^ < 0. Hence, the 
slope condition is satisfied. 

2. In the critical case d = 2/cr, the soliton power does not 
depend on /x, i.e., ^ = 0. By 14311 . the slope condition 
is violated. 

3. In the supercritical case d > 2/cr, > 0. Hence, the 
slope condition is violated. 

Thus, the slope condition is satisfied only in the subcritical 
case. 

When V = 0, the spectrum of L + is comprised of three 
(essential) parts [55], see FigureQ] 



1. A simple negative eigenvalue A m ; n < 0, with a cor- 
responding positive and radially-symmetric eigenfunc- 
tion / m i n . In 1I35IL it was shown that for power nonlin- 
earities, F(\u\) = \u\ 2a , .f m i n = u a+1 and A mm = 
a (a + 2)/i. 

2. A zero eigenvalue with multiplicity d, i.e., Ao.j = 
with eigenfunctions fj = for j — 1 , . . . , d. These 
zero eigenvalues manifest the translation invariance in 
a homogeneous medium in all d directions. 

3. A strictly positive continuous spectrum [— [i, 00). 







continuous 
spectrum^ 
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FIG. 1: The spectrum of L+ in a homogeneous medium. 

Theorem IIII.ll does not apply directly for the stability of 
solitons in homogeneous medium because Ao.j = and 
ri _ = 1. Accordingly, the notion of orbital stability must 
be modified. Indeed, by the Galilean invariance of NLS for 
Vi = V n i = 0, an arbitrarily small perturbation of a soliton 
can result in the soliton moving at small uniform speed to in- 
finity. The orbit in a homogeneous medium is thus the group 
of all translates in phase and space, i.e., {u(x — xq; /u)e 17 : 
xq e M. d , 7 € [0, 2tt) } and orbital stability is given by Def- 
inition III. II but where the infimums are taken over all 7 and 
x . 

Accordingly, Weinstein showed in ll43ll that in the case of 
homogeneous media, the spectral condition can be slightly re- 
laxed so that it is satisfied if L + has only one negative eigen- 
value and d— zero eigenvalues, associated with the transla- 
tional degrees of freedom of NLS. Hence, the spectral con- 
dition is satisfied in homogeneous media and stability is de- 
termined by the slope condition alone 14311 . In particular, 
solitons in homogeneous media with a power-law nonlinear- 
ity F(|A| 2 ) = \A\ 2(T are stable only in the subcritical case 
a < 2/d. 

B. Stability conditions in inhomogeneous media 

Below we investigate how the two stability conditions are 
affected by a potential/lattice. 

Generically, in the subcritical (d < 2/rr) and supercriti- 
cal (d > 2/cr) cases, the slope has an 0(1) magnitude in a 
homogeneous medium. Hence, a weak lattice can affect the 
magnitude of the slope but not its sign, see e.g., ll35ll . Clearly, 
a sufficiently strong lattice can alter th e sig n of the slope, see 
e.g. J36ll for the subcritical case and ll56l l57ll for the super- 
critical case. The situation is very different in the critical case 
(d = 2/cr). Indeed, since the slope is zero in a homogeneous 
medium, any potential, no matter how weak, can affect the 
sign of the slope. 

The potential can affect the spectrum of L + in two different 
ways: 1) shift the eigenvalues, and 2) open gaps (bounded- 
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intervals) in the continuous spectrum, see Figure [2] In gen- 
eral, the minimal eigenvalue of L + remains negative, i.e., 

\nin < 0' me continuous spectrum remains positive, and 
the zero eigenvalues can move either to the right or to the 
left. Hence, generically, the spectrum of L + has the following 
structure: 

1 . A simple negative eigenvalue < with a positive 
eigenfunction /^Y] > 0. 

2. Perturbed-zero eigenvalues A^ with eigenfunctions 
/j V) ,fori = l,...,d 

3. A positive continuous spectrum, sometimes with a 
band-gap structure, beginning at — fj}^ > 0. 

This structure of the spectrum was proved in [33] for solitons 
in the presence of a nonlinear lattice, i.e., Eq. (fl~|i with Vi = 0. 
For a linear lattice, the proof of the negativity of X^2 n is the 

same as in i33ll . The proof of the positivity of — /i^2 is the 
same as in H33I1 for potentials that decay to as — * oo. 



bands & 
gaps 



FIG. 2: The spectrum of L+ in an inhomogeneous medium. 



Since A^ < and the continuous spectrum is positive, 



the spectral condition (|7]l reduces to 



Ag? > 



i = i, 



(8) 



i.e., that all the perturbed-zero eigenvalues are positive. 
Generically, the equivalent spectral condition ([H} is satisfied 
when the soliton is centered at a local minimum of the poten- 
tial, but violated when the soliton is centered at a local max- 
imum or saddle point of the potential [H [H [H H3, HE 

Although generically A^ (/«) 7^ 0, there are two scenarios 
in which A V j equals zero: 

1 . The potential is invariant under a subgroup of the con- 
tinuous spatial-translation group. For example (see 
also i34lo . in a one-dimensional lattice embedded in 
2D, i.e., ^Qy'^ = 0, one has Ao,2(m) = 0- in sucn 
cases, the zero eigenvalues do not lead to instability 
for the reasons given in Section IHI Al Rather, the or- 
bit and distance function are redefined modulo the ad- 
ditional invariance, e.g., in the example above the orbit 
is {u(x, y - y ; v)e n : y G E, 7 G [0, 2tt) }. 

2. In the presence of spatial inhomogeneity, V ^ 0, Aq^ 
can cross zero as p is varied. See for example ll8~lll and 
the examples discussed in Sections rVIII Dl and llXI This 
crossing can be associated with a bifurcation and the 
existence of a new branch of solitons and an exchange 



of stability from the old to the new branch; see the sym- 
metry breaking analysis of j8lll . In such cases, stability 
and instability depend on the details of the potential and 
nonlinearity. 

In some cases, there are also positive discrete eigenvalues in 
(0, — (J,). However, these eigenvalues do not affect the orbital 
stability, since they are positive. They do play a role, however, 
in the scattering theory of solitons 1771 17811 . 

We note that in many previous studies, only the slope con- 
dition was checked for stability. As Theorem IHI. 11 shows, 
however, "ignoring" the spectral condition is justified only for 
solitons centered at lattice minima, since only then the spec- 
tral condition is satisfied. In all other cases, checking only the 
slope condition usually lead to incorrect conclusions regard- 
ing stability. 



C. Instability and collapse 

We recall that in a homogeneous medium with a power non- 
linearity, all solutions of the subcritical NLS exist globally. 
For critical and supercritical NLS there are collapsing (singu- 
lar) solutions 06111 . i.e., solutions for which j \VA(x, z)\ 2 dx 
tends to infinity in finite distance. Hence, in a homogeneous 
medium, the two phenomena of collapse and of soliton insta- 
bility appear together. In fact, the two phenomena are directly 
related, since in the critical and supercritical cases, the insta- 
bility of the solitons is manifested by the fact that they can col- 
lapse under infinitesimally small perturbations (i.e., a strong 
instability). 

As we shall see below, the situation is different in inhomo- 
geneous media. Indeed, the soliton can be unstable even if 
all solutions of the corresponding NLS exist globally. Con- 
versely, the soliton can be stable, yet undergo collapse under 
a sufficiently strong perturbation. Such results on the "de- 
coupling" of instabilityand collapse have already appeared 
in [3! [M [H [H, |6l. In all of these cases, the "decou- 
pling" is related to the absence of translation invariance. 



IV. QUALITATIVE APPROACH - CLASSIFICATION OF 
INSTABILITY DYNAMICS 



The dynamics of orbitally-stable solitons is relatively 
straightforward - the solution remains close to the unperturbed 
soliton. On the other hand, there are several possible ways for 
a soliton to become unstable: it can undergo collapse, com- 
plete diffraction, drift, breakup into separate structures, etc. 

Theorem IHI. 11 is our starting point for the classification of 
the instability dynamics, since it suggests that there are two 
independent mechanisms for (in)stability. In fact, we show 
below that the instability dynamics depends on which of the 
two conditions for stability is violated. 

As noted in Section |HI CI in a homogeneous medium with 
a power-law nonlinearity, when the slope condition is vio- 
lated, the soliton can collapse (become singular) under an 
infinitesimal perturbation. If the perturbation increases the 
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beam power, then nonlinearity dominates over diffraction so 
that the soliton amplitude becomes infinite as its width shrinks 
to zero. If the perturbation is in the "opposite direction", the 
soliton diffracts to zero, i.e., its amplitude goes to zero as its 
width becomes infinite, see e.g., Theorem 2 of B48I1 . More 
generally, in other types of nonlinearities or in the presence 
of inhomogeneities, there are cases where the slope condi- 
tion is violated but collapse is not possible (e.g., in the one- 
dimensional NLS with a saturable nonlinearity i63ll ). In such 
cases, a violation of the slope condition leads to a focusing 
instability whereby infinitesimal changes of the soliton can 
result in large changes of the beam amplitude/width, but not 
in collapse or total diffraction. Accordingly, we refer to the in- 
stability which is related to the violation of the slope condition 
as a focusing instability (rather than as a collapse instability). 

When the soliton is unstable because the spectral condition 
is violated, it undergoes a drift instability whereby infinitesi- 
mal shifts of the initial soliton location lead to a lateral move- 
ment of the soliton away from its initial location. The mathe- 
matical explanation for the drift instability is as follows. The 
spectral condition is associated with the perturbed-zero eigen- 
value Ag^ and the corresponding eigenmode fj. In the ho- 
mogeneous case, the eigenmodes fj — J^- are odd. By con- 
tinuity from the homogeneous case, the perturbed-zero eigen- 
modes / . in the presence of a potential are odd for symmet- 
ric potentials and "essentially" odd for asymmetric potentials. 
When the spectral condition is violated, these odd eigenmodes 
grow as z increases, resulting in an asymmetric distortion of 
the soliton, which gives rise to a drift of the beam away from 
its initial location. The mathematical relation between the vi- 
olation of the spectral condition and the drift instability is fur- 
ther developed in SectionlVl 

Finally, the drift dynamics also has an intuitive physical ex- 
planation. According to Fermat's Principle, light bends to- 
wards regions of higher refractive-index. Positive values of 
the potential V correspond to negative values of the refrac- 
tive index, hence, Fermat's principle implies that beams bend 
towards regions of lower potential. Moreover, since generi- 
cally, the spectral condition is satisfied for solitons centered 
at a lattice minimum but violated for solitons centered at a 
lattice maximum, one sees that the drift instability of solitons 
centered at lattice maxima and the drift stability of solitons 
centered at lattice minima is a manifestation of Fermat's prin- 
ciple. 



The quantitative relation between the value of and the 
drift rate was found analytically for the first time in l35ll for 
narrow solitons in a Kerr medium with a linear lattice. Later, 
based on the linearized NLS dynamics, it was shown in |[37ll 
that for solitons of any width, any nonlinearity and any linear 
or nonlinear potential, this quantitative relation is as follows. 
Let us define the center of mass of a perturbed soliton in the 
x ; coordinate as 



: = p / Xj\A\ 2 dx. 



(9) 



Then, by IB7I1 . the dynamics of (xj) is initially governed by 
the linear oscillator equation 

^(te> = do) 

with the initial conditions 

«j/z=o = J x 3 \A \ 2 dx/P, 



=0 = 2d ■ ImJ A*VA dx/P. 



(11) 



Here, £ ,j is the location of the lattice critical point in the jth 
direction (not to be confused with (xj) , the value of the 
center of mass at z = 0). The forcing is given by 



S~1 \(^0 (~1 

~^3 A 0,j ' u 3 



fAV) AV), 

Uoj ) Jo,j ) 
tj-iAV) f (v)y 

J0,j ) J0,j > 



(12) 



where f^y is the eigenmode of L + that corresponds to Ag^ , 
i.e., the eigenmode along the Xj direction, the operator L_ is 
given by 

L- = —A — fi — (1 — V n i(x)) F(u 2 ) + V U 

and the inner product is defined as (/, g) = J fg*dx. 

Since L_ is non-negative for positive solitons, it follows 
that Cj > 0. Therefore, when Aq^ is negative, Qj is real 

and when XqO is positive, flj is purely imaginary. Hence, by 
Eqs. ([T0l>-(112l>. it follows that the lateral dynamics of a general 
incident beam centered near a lattice minimum is 



te) x =o cos (l n jp) 



"J/*=0 



sin(|%|z), (13) 



V. QUANTITATIVE APPROACH 



As noted, the soliton is drift-unstable when XxO < but 



0,j 



drift-stable when Aq • > 0. Thus, there is a discontinuity in 

the behavior as Aq^ passes through zero. Nevertheless, one 
can expect the transition between drift instability and drift sta- 
bility to be continuous, in the sense that as Aq^ approaches 
zero from below, the rate of the drift becomes slower and 
slower. Similarly, we can expect that as Aq^- becomes more 
negative, the drift rate will increase. 



i.e., the soliton drifts along the Xj coordinate at the rate ilj. 
On the other hand, the lateral dynamics of a general incident 
beam centered near a lattice maximum is 



(xj) — (Xj) z _ cosh(Qjz) + 



dz ( X j)z=0 „; 



SI, 



sinh(Slj-z). (14) 



i.e., the soliton is pulled back towards £o,j by a restoring force 
which is proportional to fij, so that it undergoes oscillations 
around £oj in the Xj coordinate with the period |f2j|, 

As noted, the soliton is focusing-unstable when the slope 
dP/dji is non-negative, and focusing-stable when the slope 
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is negative. In a similar manner to the continuous transition 
between drift stability and instability, one can expect the tran- 
sition between focusing stability and instability to be continu- 
ous. In other words, one can expect the magnitude of the slope 
to be related to the strength of focusing stability or instability. 
At present, the quantitative relation between the magnitude of 
the slope and the strength of the stability is not known, i.e., 
we do not have a relation such as ( TTOb . However, numerical 
evidence for this link was found in several of our earlier stud- 
ies i33l [34l [35l H^]. For example, in the case of focusing - 
stable solitons that collapse under sufficiently large pertur- 
bations, it was observed that as the magnitude of the slope 
increases, the magnitude of the perturbation that is needed 
for the soliton to collapse also increases. Thus, the magni- 
tude of the slope is related to the size of the basin of stabil- 
ity l33l 1341 [3511 . In cases of focusing-stable solitons where 
collapse is not possible, when the magnitude of the slope in- 
creases, the focusing stability is stronger in the sense that for a 
given perturbation, the maximal deviation of the soliton from 
its initial amplitude decreases |36]. 



A. Physical vs. Mathematical stability 

The quantitative approach is especially important in the 
limiting cases of "weak stability/instability", i.e. when one 
is near the transition between stability and instability. For ex- 
ample, consider a soliton for which the two conditions for sta- 
bility are met, but for which or the slope are very small 
in magnitude. Such a soliton is orbitally stable, yet it can 
become unstable under perturbations which are quite small 
compared with typical perturbations that exist in experimen- 
tal setups. Hence, such a soliton is "mathematically stable" 
but "physically unstable", see e.g., [33]. Conversely, consider 
an unstable soliton for which either Aq^ is negative but very 
small in magnitude or the slope is positive but small. In this 
case, the instability develops so slowly so that it can be some- 
times neglected over the propagation distances of the experi- 
ment. Such a soliton is therefore "mathematically unstable" 
but "physically stable" JH. 



VI. GENERAL RULES 

We can summarize the results described so far by several 
general rules for stability and instability of bright positive lat- 
tice solitons. 

The qualitative approach rules are: 

QL1 Bright positive lattice solitons of NLS equations can be- 
come unstable in only two ways: focusing-instability or 
drift-instability. 

QL2 Violation of the slope condition leads to an focusing- 
instability, i.e., either initial diffraction or initial self- 
focusing. In the latter case, self-focusing can lead to 
collapse. Note, however, that for "subcritical" nonlin- 
earities, the self-focusing is arrested. 



QL3 The spectral condition is generically satisfied when the 
soliton is centered at a potential minimum and violated 
when the soliton is centered at a potential maximum or 
saddle point. 

QL4 Violation of the spectral condition leads to a drift- 
instability, i.e., an initial lateral drift of the soliton from 
the potential maximum/saddle point towards a nearby 
lattice minimum. 

The quantitative theory rules are: 

QN1 The strength of the focusing- and drift- stability and in- 
stability depends on the magnitude of the slope ^ 

and the magnitude of |Aq^ |, respectively. 

QN2 The lateral dynamics of the beam is initially given by 
Eqs. (O-flaJ. 

The above rules were previously demonstrated for ID soli- 
tons in a periodic nonlinear lattice 13311 . for an anisotropic 
2D lattice 1 34] and for several specific cases of linear lat- 
tices JlH l36ll . In this paper, we demonstrate that these rules 
apply in a general setting of dimension, nonlinearity, lin- 
ear/nonlinear lattice with any structure and for any soliton 
width. In particular, we use these general rules to explain the 
dynamics of lattice solitons in a variety of examples that were 
not studied before. 



VII. NUMERICAL METHODOLOGY 

Below we present a series of numerical computations that 
illustrate the qualitative and quantitative approaches presented 
in Sections HVllVl We present results for the 2D cubic NLS 



iA z (x, y, z) + AA + \A\ 2 A - V(x, y)A = 0, 



(15) 



with periodic lattices, lattices with a vacancy defect, and lat- 
tices with a quasicrystal structure. There are two reasons for 
the choice of the 2D cubic NLS. First, this equation enables 
us to illustrate the instability dynamics in dimensions larger 
than one, in particular, in cases where the dynamics in each 
direction is different (e.g., as for solitons centered at saddle 
points). Second, the 2D cubic NLS enables us to elucidate 
the distinction between instability and collapse. Indeed, we 
recall that a necessary condition for collapse in the 2D cubic 
NLS is that the power of the beam exceeds the critical power 

P c « 11.7 & 

We first compute the soliton profile by solving Eq. (f2| us- 
ing the spectral renormalization method 16411 . Once the soli- 
tons are computed for a range of values of p, the slope condi- 
tion © is straightforward to check. In order to check the spec- 
tral condition (0, the perturbed-zero eigenvalues A^ (and 
the corresponding eigenfunctions fj) of the discrete approxi- 
mation of the operator L + are computed using the numerical 
method presented in [35, Appendix D]. The value of flj is 
calculated from Eq. ( fT2l by inversion of the discrete approxi- 
mation of the operator L_. 
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Eq. ( TT3T > is solved using an explicit Runge-Kutta four-order 
finite-difference scheme. Following ll33i[34l[35i[36ll . the ini- 
tial conditions are taken to be the unperturbed lattice soliton 
u(x, y) with either 

1. a small power perturbation, i.e., 

Aq(x, y) = Vl + cu(x, y), (16) 



(a) (b) 

' ;•»;♦; " vv in 

y 0* l^^r- saddle 

!♦!•!♦!♦! • " 

-1 x 1 -0 5 -0.5 a 



where c is a small constant that expresses the excess 
power of the input beam above that of the unperturbed 
soliton, or 

2. a small lateral shift, i.e., 



FIG. 3: (Color online) The sinusoidal square lattice given by Eq. l |19t 
with Vb = 5. (a) Top view, (b) Side view. The solitons investigated 
below are centered at the lattice maximum (0,0), lattice minimum 
(0.25,0.25), and saddle point (0.25,0). 



Ao(x,y) = u(x - Ax ,y~ Ay ), 



(17) 



A. Solitons at lattice minima 



where Axo and Ayo are small compared with the char- 
acteristic length-scale (e.g., period) of the potential. 

The motivation for this choice of perturbations is that each 
perturbation predominantly excites only one type of instabil- 
ity. Indeed, by Eq. (fT0b-(fTTT>. it is easy to verify that under 
a power perturbation ([ToT i the center of mass will remain at 
its initial location ("cf. ll33l l34l l36To . i.e., no lateral drift will 
occur. In this case, only an focusing instability is possible. 
On the other hand, the asymmetric perturbation dTvT > will pre- 
dominantly excite a drift instability (but if the soliton is drift- 
stable, this perturbation can excite an focusing instability, see 
Figure |6]l. 

The advantage of the perturbations (TT&b - (TTTb over adding 
random noise to the input soliton is that they allow us to con- 
trol the type of instability that is excited. Moreover, grid con- 
vergence tests are also simpler. Once the NLS solution is com- 
puted, it is checked for focusing and drift instabilities by mon- 
itoring the evolution of the normalized peak intensity 



I(z) := 



max x , y \A(x, y, 



\A (x,y)\ 2 ' 
and of the center of mass (|9), respectively. 

VHL PERIODIC SQUARE LATTICES 

We first choose the sinusoidal square lattice 



(18) 



V(x,y) 



To 



cos 2 (27ra;) 



'■(2ny)] 



(19) 



which is depicted in Figure [3] We consider this to be the sim- 
plest 2D periodic potential, as all the local extrema are also 
global extrema. This lattice can be created through interfer- 
ence of two pairs of counter-propagating plane waves, and is 
standard in experimental setups, see, e.g., |j65ll66ll . The stabil- 
ity and instability dynamics are investigated below for solitons 
centered at the lattice maxima, minima, and saddle points, see 
Figure EJb). 



We first investigate solitons centered at the lattice minimum 
(xo,yo) = (0.25,0.25). Figure|4fa) shows that the power of 
solitons at lattice minima is below the critical power for col- 
lapse, i.e., P(p) < P c ~ 11-7 for all fi. As the soliton be- 
comes narrower (/i — > — oo), the soliton power approaches P c 
from below (as was shown numerically in Il40ll for this lattice 
and analytically in ll3~5ll for any linear lattice). In addition, as 
the soliton becomes wider ([i — » /iS^, the edge of the first 
band), its power approaches P. from below (rather than be- 
comes infinite, as implied in IHOIl ). see also 19111 . The minimal 
power is obtained at p = /i m = —10. The power curve thus 
has a stable branch for narrow solitons (— oo < /i < p m ) 
where the slope condition is satisfied, and an unstable branch 
for wide solitons (p m < p < /U^g) where the slope condi- 
tion is violated. Therefore, wide solitons should be focusing- 
unstable while narrow solitons should be focusing-stable. Fig- 
ure[4|b) shows that, as expected for solitons at lattice minima, 
Ag 1 "* = A Q 2 ^ > for all ji. Hence, the spectral condition is 
fulfilled. Consequently, solitons at lattice minima should not 
experience a drift instability. 

In order to excite the focusing instability alone, we add 
to the soliton a small power perturbation, see Eq. ( fT6l) . We 
contrast the dynamics in a neighborhood of stable and unsta- 
ble solitons by choosing two solitons with the same power 
(P = 0.98P C ), from the stable branch (p = —31) and from 
the unstable branch (p = —3). We perturb these solitons with 
the same power perturbations (c = 0.5%, 1%, 2%). 

When c = 0.5% and 1%, the input power is below the 
threshold for collapse (P < P c ). In these cases, the self- 
focusing process is arrested and, during further propagation, 
the normalized peak intensity undergoes oscillations (see Fig- 
ures |3 a) and (b)). For a given perturbation, the oscillations 
are significantly smaller for the stable soliton compared with 
the unstable soliton. 

When c = 2.5%, the input power is above the threshold for 
collapse (P > P c ) and the solutions undergo collapse. There- 
fore, for such large perturbations, collapse occurs for both sta- 
ble and unstable solitons, i.e., even when both the slope and 
spectral conditions are fulfilled. This shows yet again that in 
an inhomogeneous medium, collapse and instability are not 
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necessarily correlated. 
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FIG. 4: (Color online) (a) Power, and (b) perturbed-zero eigenval- 
ues, as functions of the propagation constant, for solitons centered 
at a maximum (blue, dashes) and minimum (red, dots) of the lat- 
tice U9\ with Vo — 5. Also shown are the corresponding lines for 
the homogeneous NLS equation (solid, green). The circles (black) 
correspond to the values used in Figs.| 
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FIG. 5: (Color online) Normalized peak intensity d!8t of solutions 
of Eq. J 15b with the periodic lattice J19l > with Vq = 5. Initial condi- 
tions are power-perturbed solitons [see Eq. d!6b l centered at a lattice 
minimum: (a) Soliton from the stable branch (/i = —31); (b) Soliton 
from the unstable branch {p — —3). Input powers are 0.5% (red 
dots), 1% (blue dashes), and 2.5% (solid green) above the soliton 
power. 



instability is excited and the soliton amplitude decreases (as 
its width increases), see Figure |6fb2). Obviously, once the 
soliton amplitude changes significantly, the theoretical pre- 
diction for the lateral dynamics is no longer valid. In order 
to be convinced that the initial instability in this case is of an 
focusing-type rather than drift-type, we note that for fi = —31 
for which the slope condition is satisfied, the soliton remains 
focusing-stable, see Figure[6];a2). 
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FIG. 6: (Color online) Dynamics of solutions of Eq. d!5b with the 
periodic lattice d!9b with Vo = 5. Initial conditions are position- 
shifted solitons [see Eq. d!7b l centered at a lattice minimum, with 
(Axo, Ayo) = (0, 0.04). (al) Center of mass in the y coordinate 
(blue, solid line) and analytical prediction [Eq. J20t . red, dashes] for 
fi = —31; (a2) Normalized peak intensity d!8b for /i = —31; (bl) 
and (b2) are the same as (al) and (a2), but for fi — —3. 



In order to confirm that solitons centered at a lattice min- 
imum do not undergo a drift instability, we shift the soli- 
ton slightly upward by using the initial condition (Fi\ with 
(Axo,Ayo) = (0,0.04). Under this perturbation, the solu- 
tion of Eq. (fTOb is 

(x}=0, (y) = Ay -cos(\n y \z). (20) 

In addition, by Eq. ( fT2l . Sl y w 11.12i for /.i = —31 and VL y w 
2.58i for /i = —3. Figure|6|al) shows that for fi = —31, the 
center of mass in the y-direction of the position-shifted soliton 
follows the theoretical prediction d20l ) accurately over several 
oscillations. In addition, the center of mass in the x-direction 
remain at x = (data not shown), in agreement with Eq. d20l ). 
Thus, the soliton is indeed drift-stable. 

The situation is more complex for /i — — 3. In this case, the 
position-shifted soliton follows the theoretical prediction ( f20b 
over more than 2 diffraction lengths (i.e., for z > zq where 
zq 1), but then deviates from it, see Figure |6|bl). The 
reason for this instability is that for fi = —3, the slope con- 
dition is violated. Since the position-shifted initial condition 
can also be viewed as an asymmetric amplitude power pertur- 
bation AA = u(x — Axq, y — Ayo) — u(x, y), an focusing 



B. Solitons at lattice maxima 

We now investigate solitons centered at the lattice maxi- 
mum (xo,yo) = (0,0). Figure [4] shows that in general, soli- 
tons at lattice maxima have the opposite stability characteris- 
tics compared with those of solitons centered at lattice min- 
ima: The slope condition is violated for narrow solitons and 
satisfied for wide solitons, the power is above P c ll92ll . and the 
perturbed-zero eigenvalues Ag? are always negative. Inter- 
estingly, for the specific choice of the lattice ( fl9] >, the powers 
and perturbed-zero eigenvalues at lattice maxima and minima 
are approximately, but not exactly, images of each other with 
respect to the case of a homogeneous medium. 

The negativity of the perturbed-zero eigenvalues implies 
that solitons centered at a lattice maximum undergo a drift 
instability (see Figure |HJb)). However, if the initial condi- 
tion is subject to a power perturbation, see Eq. (fToT i. then no 
drift occurs. In this case, stability is determined by the slope 
condition. For example, Figure Q shows the dynamics of a 
power-perturbed wide soliton for which the slope condition 
is satisfied. When the soliton's input power is increased by 
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0.5%, the solution undergoes small focusing-defocusing os- 
cillations, as in Figure |5f a), i.e., it is stable under symmetric 
perturbations. When the soliton's input power is increased by 
1%, the perturbation exceeds the "basin of stability" of the 
soliton 1 35] and the soliton undergoes collapse. These results 
again demonstrate that collapse and instability are indepen- 
dent phenomena. 

If the initial condition is asymmetric with respect to the lat- 
tice maximum, the soliton will undergo a drift instability. In 
Figure [8] we excite this instability with a small upward shift, 
namely, Eq. ^ with (Ax ,Ay ) = (0,0.02). Under this 
perturbation, the solution of Eq. dTob is 



(a) 



(b) 



(y) = Ayo ■ cosh(n y z). 



(21) 



with Q y ps 3.9. In the initial stage of the propagation 
(z < 0.5) the soliton drifts toward the lattice minimum - 
precisely following the asymptotic prediction ( 1201 ). see Fig- 
ure [Sf a), but the soliton's amplitude is almost constant), see 
Figure [8fb). During the second stage of the propagation 
(0.5 < z < 0.99) the soliton drifts somewhat beyond the 
lattice minimum as it begins to undergo self-focusing. In the 
final stage (0.99 < z < 1) the soliton undergoes collapse (Fig- 
ure[8|b)). The global dynamics can be understood in terms of 
the stability conditions for solitons centered at lattice minima 
and maxima as follows. The initial soliton, which is centered 
at a lattice maximum, satisfies the slope condition but violates 
the spectral condition. Consistent with these traits, the soli- 
ton is focusing-stable but undergoes a drift instability. As the 
soliton gets closer to the lattice minimum, it can be viewed 
as a perturbed soliton centered at the lattice minimum, for 
which the spectral condition is fulfilled and the soliton power 
is below P c (see Figure |4|b)). Indeed, at this stage, the drift 
is arrested because the beam is being attracted back towards 
the lattice minimum. Moreover, the beam now is a strongly 
power-perturbed soliton, since the beam power (ps 1.03P c ) is 
w 6% above the power of the soliton at a lattice minimum. 
Hence, in a similar manner to the results of Figure |5ja), the 
perturbation exceeds the "basin of stability" and the soliton 
undergoes collapse. 
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FIG. 7: (Color online) Same as Figure^a) for a soliton at a lattice 
maximum with fi = —5 (stable branch) and input power that is 0.5% 
(red dots) and 1% (blue dashes) above the soliton power. 
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FIG. 8: (Color online) Dynamics of a soliton at a lattice maxi- 
mum with /j, = — 5, which is position-shifted according to d!7l l with 
(Axo, A?/o) ~ (0,0.02). (a) Center of mass in the y coordinate 
(blue, dashes) and the analytical prediction (Eq. ill I with Q y w 3.9, 
solid black). Location of lattice minimum and maxima are denoted 
by thin magenta and black horizontal lines, respectively, (b) Normal- 
ized peak intensity. 



C. Solitons at a saddle point 

From the didactic point of view, it is interesting also to con- 
sider solitons centered at a saddle point since they exhibit a 
combination of the features of solitons at lattice minima and 
maxima. To show this, we compute solitons centered at the 
saddle point (xo, ya) = (0.25, 0) of the lattice ( fT9l . 

Figure |9[a) shows that the zero eigenvalues bifurcate into 
Ag 1 '' > on the stable x-direction, i.e., along direction in 

(2) 

which the saddle is a minimum, and to Ag < on the unsta- 
ble y-direction, where the saddle in a maximum. 

The opposite signs of the perturbed-zero eigenvalues im- 
ply a different dynamics in each of these directions. In or- 
der to excite only the drift instability, we solve Eq. (TT3T > with 
fj, = —12 which belongs to the focusing-stable branch (see 
Figure|9jb2)). For this value of fi, the perturbed zero eigenval- 



ues are A^ = 1.7 and A Q 2J = -1.8. By <fl2i the theoretical 
prediction for the oscillation period is Q x = \7i\ = 7 whereas 
the drift rate is Sl y = 7.2. Hence, the theoretical prediction 
for the dynamics of the center of mass is 

(x) fa 0.25 + Ax ■ cos(7z), (y) ss Ay • cosh(7.2z). 

Indeed, a shift in the x direction (Axq, Ayo) ~ (0.0156, 0) 
leads to oscillation in the x-direction (Figure [lOfa)) while 
(y) (Figure Q/Jtb)) and the amplitude (Figure [TOfc)) are un- 
changed. On the other hand, a shift in the y direction 
(Axo, Ayo) w (0,0.0156) leads to a drift instability in the 
^-direction (Figure [TOl'b)) but has no effect on (x) (Fig- 
ure [Tol a)). In both the stable and unstable directions, the cen- 
ter of mass follows the analytical prediction remarkably well. 
Figure [TOl c) also shows that once the soliton drifts beyond 
the lattice minimum, the beam undergoes collapse. This can 
be understood using the same reasoning used for solitons that 
drift from a lattice maximum (see explanation for Figure [8] in 
Section lVrrTBl . 

We also note that for the specific choice of the lattice (\9l , 
the values of the perturbed-zero eigenvalues in the stable and 
unstable directions are nearly indistinguishable from those 



,(2) 



of the perturbed-zero eigenvalues that correspond to solitons 
centered at a lattice minimum and maximum, respectively. 
This can be understood by rewriting the lattice ( fT9l l as 

V(x, y) = -J [1 - cos 2 (2tt(2; - 0.25)) + cos 2 (27ry)] . 

(22) 

Thus, apart from the constant part (i.e., the first term), the 
difference between the lattices is the sign before the x- 
component of the lattice. In that sense, in the x direction, 
the saddle point is equivalent to a maximum point, hence, the 
similarity between the eigenvalues. Another consequence of 
the x — y symmetry of the lattice ( fl"9l > is that the soliton has 
approximately the critical power P c for all fi, i.e., P(fi) ~ P c , 
which is approximately the average of the powers of solitons 
at maxima and minima (see Figure |9jb 1 ,b2)). As noted be- 
fore, this will be no longer true if the lattice changes in the x 
and y directions will no longer be equal. 
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FIG. 9: (Color online) (a) The perturbed-zero eigenvalues at the sad- 
dle point. One eigenvalue is shifted to positive values (magenta), and 
is indistinguishable from the eigenvalue at lattice minima (red); one 
eigenvalue is shifted to negative values (black), and is indistinguish- 
able from the eigenvalue at lattice maxima (blue), (bl) Same data as 
in Figure|4la), with the addition of data for solitons centered at a sad- 
dle point of the lattice (black, dash-dots). (b2) same as (bl) showing 
only the data for solitons centered at a saddle point (black, dash-dots) 
and for the homogeneous medium soliton (greed line). 



(a) (b) (c) 




FIG. 10: (Color online) Dynamics of a soliton centered at a saddle 
point (x ,y ) = (0.25,0) of the lattice {[9} with y, = -12 and 
shifts along: (/) the stable x direction [(Aa;o, Aj/o) ~ (0.0156, 0), 
red dots], (ii) the unstable y direction [(Aso, Ayo) ~ (0, 0.0156), 
green dashes], (in) the diagonal direction [(Axo,Ayo) ~ 
(0.0156, 0.0156), blue, dash-dots], (a) Center of mass (x). (b) Cen- 
ter of mass (y). (c) Normalized peak intensity. 

If we apply perturbations in the stable and unstable direc- 
tions simultaneously (Axo, Ayo) ~ (0.0156, 0.0156), the dy- 
namics in each coordinate is nearly identical to the dynamics 
when the perturbation was applied just in that direction. Thus, 
there is a "decoupling" between the (lateral dynamics in the) 
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x and y directions. Indeed, this decoupling follows directly 
fromEq. ((101). 

D. Solitons at a shallow-maximum 

We now consider solitons of the periodic potential 

V{x, y) = [2 cos(27ra) + 2 cos(27ry) + l] 2 , (23) 
25 

where Vq = 5 and the normalization by 25 implies that Vq = 
max Ij9 V(x, y). Unlike the lattice ( fl"9l >, the lattice d23l also 
has shallow local maxima that are not global maxima [e.g., at 
(0.5,0.5)]. 

The stability and instability dynamics of solitons centered 
at global minima, maxima and saddle points of the lattice are 
similar to the case of the lattice (19\ , which was already stud- 
ied. Hence, we focus only on the stability of solitons centered 
at a shallow maximum. 

Since the lattice is invariant under a 90° rotation, the 
perturbed-zero eigenvalues are equal, i.e., Aq = A 2 \ How- 
ever, unlike solitons centered at a global maximum, the cor- 
responding perturbed-zero eigenvalues are negative only for 
very negative values of fi (narrow beams) but become positive 
for values of n near the band edge \ibe (wide beams), see Fig- 
ure [T2T b). The reason for the positivity of Aq = Aq 2 ' despite 
being centered at a lattice maximum is as follows. For narrow 
solitons, the region where the "bulk of the beam" is located is 
of higher values of the potential compared with the immedi- 
ate surrounding, hence, the solitons "feel" an effective lattice 
maximum. On the other hand, for wider solitons, the "bulk of 
the beam" is centered mostly at the shallow lattice maximum 
and the surrounding lower potential regions. Hence, although 
the very center of the soliton is at the shallow lattice maxi- 
mum, these solitons are effectively centered at the lattice min- 
imum with respect to the nearest global lattice maxima (see 
also HH, Section 4.5). The transition of the qualitative sta- 
bility properties between narrow and wide solitons described 
above occurs when the soliton's width is on the order of the 
lattice period. As noted in Section IIII Bl the stability at the 
transition points where Aq^ = or ^ =0 requires a specific 
study. Similarly, a comparison of Figure [T27 a) and FigureUa) 
shows that the P(/i) reflects the transition between properties 
which are characteristic to solitons centered at lattice maxima 
and minima. Indeed, for narrow solitons (/i — * — oo) is similar 
to the power of solitons centered at a global maximum, i.e., the 
power is above critical and the slope is positive. On the other 
hand, P{p) curve for wide solitons (fi — ► /i^2) is similar to 
the power of solitons centered at a (simple) lattice minimum, 
i.e., the power is below critical and the slope is positive too. 

Numerical simulations (Figure Qj) demonstrate this transi- 
tion. For a narrow soliton (fi = —12), the theoretical predic- 
tion for the dynamics of the center of mass is (x) = 0.5 + 
Ax cosh(4.14z) and (y) = 0.5 + Ay cosh(4.14z). Indeed, 
the narrow soliton drifts away from the shallow maximum to- 
ward the nearby (global) lattice minimum (Figure [T3l al)) and 
then undergoes collapse (Figure [T3l' a2)). This dynamics is 
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similar to that of solitons centered near lattice maximum or 
a saddle of a the lattice (O, see Sections [VIII B| and lVIII CI 
On the other hand, for the wide soliton (/i = —2), the the- 
oretical prediction for the dynamics of the center of mass is 
(x) = 0.5 + Axo cos(1.6z) and (y) = 0.5 + Ay cos(1.6z). 
Indeed, this soliton remains stable, undergoing small position 
oscillations around the shallow maximum (FigureQjJb)). This 
dynamics is the same as for solitons centered at a minimum of 
the lattice ( fT9] l, see Figure[6fa). As in previous examples, the 
numerical results are in excellent agreement with the analytic 
prediction dT0b-dT2l). 



(b) 
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IX. PERIODIC LATTICES WITH DEFECTS 

Defects play a very important role in energy propagation 
through inhomogeneous structures. They arise due to imper- 
fections in natural or fabricated media. They are also often 
specifically designed to influence the propagation. 

Solitons in periodic lattices with defects have drawn much 
attention both experimentally and theoretically; see, for ex- 
ample, lfl6ll67ll68l l79tl . The complexity of the lattice details 
offers an opportunity to demonstrate the relative ease of apply- 
ing the stability/dynamics criteria to predict and decipher the 
soliton dynamics in them. As an example, we study lattices 
with a point defect. Our analysis can also extend to different 
types of defects such as line defects, see e.g. |lql . 

We consider the lattice ( T231 



x - y 



V(x, y) = — 2 cos(27ra) + 2 cos(27n/) + e l0{x > y] 
25 

where the phase function 8(x, y) is given by 



, (24) 



FIG. 1 1 : (Color online) The shallow maximum periodic lattice given 
by Eq. J23t with Vb = 5 . (a) Top view, (b) Side view, (c) Cross 
section along the line x = y. 




FIG. 12: (Color online) Same as Figure [4] for solitons centered 
at a shallow local maximum of the shallow-maximum periodic lat- 
tice <[23T >- 
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FIG. 13: (Color online) Dynamics of a perturbed soliton at shallow- 
maximum periodic lattice d23t with a narrow soliton [(al) and (a2) 
with /i = —12] and a wide soliton [(b) with /i = —2], and using 
(Ax ,Ay ) = (0.05,0.05). (al) Center of mass (x) = (y) of the 
narrow soliton (blue, dashes) and the analytical prediction (red dots). 
(a2) Normalized peak intensity of the narrow soliton. (b) Same as 
(al) for the wide soliton. 
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see Figure [141 and also fltill . Compared with the shallow- 
maximum periodic lattice ( 1231 ). here the constant (DC) com- 
ponent (the third term in the lattice) attains a phase distortion 
which creates an (effective) vacancy defect at (0, 0), which is 
a shallow-maximum. Further, far away from the origin, the 
potential (1241 is locally similar to the shallow-maximum peri- 
odic lattice ( 1231 ). This is a generic example of a point defect, 
as opposed to a line defect [69]. In what follows, we consider 
solitons centered at the vacancy defect (x , y ) = (0,0). 

The stability properties of solitons in the shallow-maximum 
periodic ( 1231 ) and vacancy-defect d2"4"i i lattices are strikingly 
similar, as can be seen from Figs. Q~2] and Q3] In both cases, 
there is a marked transition between narrow and wide solitons 
and this transition occurs when the soliton width is of the or- 
der of the lattice period. Indeed, numerical simulations show 
that the dynamics of perturbed solitons is qualitatively similar 
in both cases - compare Figures [T3l and [T6l We do note that 
unlike the shallow-maximum periodic lattice, the perturbed- 
zero eigenvalues of the vacancy lattice bifurcate into differ- 
ent, though similar, values. The reason for this is the phase 
function d25l ) is not invariant by 90° rotations. 

Inspecting the lattice surfaces (Figures [TT] and H4l . it is 
clearly seen that the reason for the similarity between the 
shallow-maximum periodic and vacancy lattices is that the va- 
cant site is essentially a shallow local maximum itself - and 
only a bit shallower than those of the shallow-maximum peri- 
odic lattice (see Figure[T4"l). 

In Figure [17] we give a detailed graphical illustration of a 
typical instability dynamics due to a violation of the spectral 
condition. Figure [PTl aHc) show contours of the soliton pro- 
files superposed on the contour plot of the lattice. It can be 
seen that as a result of the initial position shift, the soliton 
drifts towards the lattice minimum and that it self-focuses at 
the same time. Figure \FJ\ d) shows the trajectory of the beam 
across the lattice. In addition, Figure[T7]'e) shows the center of 



12 



mass dynamics as a function of the intensity I(z). This shows 
that initially, the perturbed soliton undergoes a drift instability 
with little self-focusing, but that once the collapse accelerates, 
it is so fast so that the drift dynamics becomes negligible. 





FIG. 14: (Color online) Same as Figure QT] for solitons centered at 
the "vacancy" of the lattice ( 124b . 
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(Color online) (a)-(c): contours of the inten- 
|2 (blue) superimposed on the vacancy lattice (green) 



FIG. 17: 

sity \u(x,y,z) 

with initial conditions corresponding to the mode with n — — 8 that 
is initially shifted in the (a;, y) plane to (Also, Ayo) = (0.05, 0.1), 
i.e., at an angle of 63° to the y axis, (a) z = 0, 7 « 1, (b) 
z = 0.51, 7 w 2.18, (c) 2 = 0.63, 7 « 11.1. (d) Center of 
mass dynamics (black curve) and the analytical prediction (magenta, 
dashes) superimposed on the contours of the potential (green), (e) 
(x) (blue, solid) and (y) (red, dashes) as functions of I(z). Circles 
(black) correspond to the z-slices shown in (a)-(c). 



FIG. 15: (Color online) Same as Figure|4]for solitons at the vacancy 
of the lattice ( 124b . (b) The perturbed-zero eigenvalues Aq 1 ''' 2 ' are 
slightly different from each other. The circles (black) correspond to 
the values used in Figure [16] 
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FIG. 16: (Color online) Same as Figure[73lbut for the vacancy lat- 
tice ((24). Here fl x » 3 in (a2) and Q x w 1.09i in (b). In both cases 
the (y) dynamics (not shown) is similar (but not identical) to the (x) 
dynamics. 



X. QUASICRYSTAL LATTICES 

Next, we investigate solitons in quasicrystal lattices. Such 
lattices appear naturally in certain molecules | [ 7 0| , [7 111. have 
been investigated in optics lfl6l l23l l24l l25l l26ll and in 
BEC B72I1 . and can be formed optically by the far-field diffrac- 
tion pattern of a mask with point-apertures that are located on 
the N vertices of a regula r polyg on, or equivalently, by the 



sum of N plane waves (cf. Ilq,l73ll ) with wavevectors (k x , k y ) 



whose directions are equally distributed over the unit circle. 
The corresponding potential is given by 



V{x,y) 



where (k ( x n) , k ( y n) ) = (K cx>s(2nn/N),K8ia(2im/N)) (93|]. 
The normalization by iV 2 implies that Vb = max x ,y V(x, y). 
The potential (f26b with N = 2, 3, 4, 6 yields periodic lattices. 
All other values of N correspond to quasicrystals, which have 
a local symmetry around the origin and long-range order, but, 
unlike periodic crystals, are not invariant under spatial trans- 
lation ieI. 

We first consider the case N = 5 (a 5-fold symmetric "Pen- 
rose" quasicrystal) for solitons centered at the lattice maxi- 
mum (xQ,yo) — (0, 0), see FigurefTSI Since the soliton profile 
and stability are affected mostly by the lattice landscape near 
its center, we can expect the stability properties of the Pen- 
rose lattice soliton at (0, 0) to be qualitatively the same as for 
a soliton at a lattice maximum of a periodic lattice. Indeed, 
Figure [19] reveals the typical stability properties of solitons 
centered at a lattice maximum: An focusing-unstable branch 
for narrow solitons, an focusing-stable branch for wider soli- 
tons and negative perturbed zero-eigenvalues (compare e.g. 
with Figure [8]). Therefore, the Penrose soliton will drift from 
the lattice maximum under asymmetric perturbations and if 
the soliton is sufficiently narrow, it can also undergo collapse. 

Figure[l9]presents also the data for a perfectly periodic lat- 
tice (N = 4) and for a higher-order quasicrystal (N = 11). 



N 2 



N-l 

e l( > k '- l}x+kl y n}y) 

71 = 



(26) 
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One can see that the stability properties in these lattices is 
qualitatively similar to the N — 5 case. The only marked 
difference as N increases is that the soliton's power becomes 
larger for a given 

These results show that in contrast to the significant effect 
of the quasi-periodicity on the dynamics of linear waves (com- 
pared with the effect of perfect periodicity l24ll ). the effect of 
quasi-periodicity on the dynamics of solitons is small. 



(a) (b) 




FIG. 18: (Color online) Same as Figure [TTT a)+(b) for the Penrose 
quasicrystal lattice given by Eq. \2b\ with N = 5 and Vo = 5. 
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FIG. 19: (Color online) Same as Figure|4]for solitons at the maxima 
of the lattices i26i with N = 4 (periodic lattice, dashed blue line), 
TV = 5 (Penrose quasicrystal lattice, dash-dotted red line), N = 11 
(higher-order quasicrystal lattice, dotted black line), and the homo- 
geneous NLS soliton (solid green line). 



XI. SINGLE WAVEGUIDE POTENTIALS 

So far we studied periodic, periodic potentials with defects 
and quasiperiodic potentials. However, our theory can be ap- 
plied to other types of potentials. Indeed, let us consider lo- 
calized potentials, such as single or multiple waveguide po- 
tentials, for which the potential decays to zero at infinity. For 
such potentials, there are two limits of interest. The first limit 
is of solitons which are much wider than the width of the po- 
tential. In this case, the potential can be approximated as a 
point defect in an homogeneous medium. Then, the dynamics 
is governed by 

iA z (x, z) + AA + \A\ 2a A - jS{x)A = 0, (27) 

where 7 is a real constant. In [36], the qualitative and quan- 
titative stability approaches were applied to Eq. (|27T i in one 
transverse dimension. 



The second limit is of solitons which are much narrower 
than the width of the potential. In this case, only the local 
variation of the potential affects the soliton profile and stabil- 
ity. Hence, the potential can be expanded as 

V(x) = V{0) + ^V"(0)x 2 + ■■■ . 

The qualitative and quantitative stability approaches were ap- 
plied to this case in 113511 . 

In 1I35I HfJ, the profiles, power slope and perturbed-zero 
eigenvalues were computed analytically (exactly or asymptot- 
ically). It was proved that the perturbed-zero eigenvalues are 
negative for solitons centered at lattice maxima (repulsive po- 
tential) and are positive for solitons centered at lattice minima 
(attractive potential). Hence, in the latter case, stability is de- 
termined by the slope condition. In those two studies, detailed 
numerical simulations confirmed the validity of the qualitative 
and quantitative approaches. Hence, we do not present a sys- 
tematic stability study for localized potentials. 



XII. FINAL REMARKS 

In this paper, we presented a unified approach for analyz- 
ing the stability and instability dynamics of positive bright 
solitons. This approach consists of a qualitative character- 
ization of the type of instability, and a quantitative estima- 
tion of the instability growth rate and the strength of stability. 
This approach was summarized by several rules (Section [VTb 
and applied to a variety of numerical examples (Sections lVIIIf 
IXl ). thus revealing the similarity between a variety of physical 
configurations which, a priori, look very different from each 
other. In that sense, our approach differs from most previous 
studies which considered a specific physical configuration. 

One aspect which was emphasized in the numerical exam- 
ples is the excellent agreement between direct numerical sim- 
ulations of the NLS and the reduced equations for the center 
of mass (lateral) dynamics, Eqs. (TT0b-(fT2b. Different reduced 
equations for the lateral dynamics were previously derived un- 
der the assumption that the beam remains close to the initial 
soliton profile (see e.g. l75ll ) or by allowing the soliton param- 
eters to evolve with propagation distance (see e.g., 117611 and 
references therein). These approaches, as well as ours, are 
valid only as long as the beam profile remains close to a soli- 
ton profile. However, unlike previous approaches, Eqs. ( fTOb - 
(fT2b incorporate linear stability (spectral) information into the 
center of mass dynamics. Thus our approach shows that the 
beam profile evolves as a soliton perturbed by the eigenfunc- 
tion . The validity of this perturbation analysis is evident 
from the excellent comparison between the reduced Eqs. ( fTOb - 
( fT2b and numerical simulations for a variety of lattice types. 
To the best of our knowledge, such an agreement was not 
achieved with the previous approaches. 

The numerical examples in this paper were for two- 
dimensional Kerr media with various linear lattices. Together 
with our previous studies which were done for narrow solitons 
in any dimension l35ll . a linear delta-function potential j36Tl 
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and for nonlinear lattices l33i l34ll . there is a strong numeri- 
cal evidence that our qualitative and quantitative approaches 
apply to positive solitons in any dimension, any type of non- 
linearity of type F(|j4| 2 ) (e.g., saturable) as well as for other 
lattice configurations, e.g., "surface" or "corner" solitons IU2I1 . 

Theorem lHI.ll as well as the qualitative and quantitative ap- 
proaches apply also for the e?-dimensional discrete NLS. This 
equation is obtained from Eq. ( TT3T > by replacing A by the dif- 
ference Laplacian operator on a discrete lattice and V by a 
potential defined at discrete lattice sites. This model was ex- 
tensively studied, mostly for periodic lattices, see e.g., for ID 
and 2D discrete NLS equation with cubic nonlinearity (see 
e.g., JO, EH I3S IH), saturable nonlinearity (see e.g., H), 
cubic-quintic nonlinearity (see e.g., l85lo . General results 
on existence and stability of solitons in d— dimensions with 
power nonlinearities appear in 11821 18611 . Indeed, for the dis- 
crete NLS, the operator L + does not generically have a zero 
eigenvalue due to absence of continuous translation symme- 
try, and the continuous spectrum is a bounded interval, start- 



ing at the soliton frequency, —\i [86]. However, these changes 
in the spectrum do not affect the stability theory, the possible 
types of instabilities and the analysis of their strength. 

As noted, our analysis shows that for positive bright soli- 
tons, only two types of instabilities are possible - focusing 
instability or drift instability. Other types of instabilities may 
appear, but only for non-positive solitons (e.g., gap solitons or 
vortex solitons). A formulation of a qualitative and quantita- 
tive theories for such solitons requires further study. 
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